A microscopic approach to phase transitions in quantum systems 
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We present a new theoretical approach for the study of the phase diagram of interacting quan- 
tum particles: bosons, fermions or spins. In the neighborhood of a phase transition, the expected 
renormalization group structure is recovered both near the upper and lower critical dimension. In- 
formation on the microscopic hamiltonian is also retained and no mapping to effective field theories 
is needed. A simple approximation to our formally exact equations is studied for the spin-5 1 Heisen- 
, berg model in three dimensions where explicit results for critical exponents, critical temperature 

D ■ and coexistence curve are obtained. 

^: 

i Several physical systems, ranging from magnets to superfluids and superconductors, display rich phase diagrams 

in a temperature regime where quantum effects cannot be neglected. Different scenarios, characterized by competing 
order parameters and zero temperature phase transitions have been recently advocated also in the framework of high 
temperature superconductivity where antiferromagnetic order, Cooper pairing and, possibly, phase separation are 
at play in the same region of the phase diagram Q. A satisfactory understanding of phase transitions in quantum 
models has been attained years ago through the seminal work by Hertz Q who showed that, at low energy and long 
wavelengths, quantum models may be described by a suitable classical action. However, a quantitative theory of the 
thermodynamic behavior is still lacking and we mostly rely on mean field approaches or weak coupling renormalization 
group (RG) calculations || , applied to quantum systems via the mapping to the appropriate effective field theory. In 
particular, the interplay between thermal and quantum fluctuations is expected to give rise to crossover phenomena 
whose extent strongly depends on the microscopic features of the system. Even for the most extensively studied 
models, like the Heisenberg antiferromagnet, our knowledge of the phase diagram is in fact limited, and the first 
precise finite temperature simulation attempting to fill this gap has become available only recently Q. By contrast, 
in classical models, numerical simulations are quite efficient even in the neighborhood of critical points [^) and, from 
(Zf-j . the analytical side, microscopic approaches especially devised for the quantitative description of the phase diagram 
£C) ' of classical fluids and magnets are available. For instance, the hierarchical reference theory of fluids (HRT) Q has 
proven quite accurate in locating the phase transition lines both in lattice and in continuous models. 

In this Letter we sketch the derivation of the quantum hierarchical reference theory of fluids (QHRT) which we then 
apply to the Heisenberg antiferromagnet. We will demonstrate that the known renormalization group equations near 
four and near two dimensions are naturally recovered within our approach, which therefore unifies two complimentary 
. techniques. On approaching the critical point, the spin velocity vanishes according to the expected dynamical critical 
exponent for an antiferromagnet. Finally, the phase diagram of this model in three dimensions is computed by 
numerical integration of a simple approximation to the the QHRT equations, providing a concrete application of our 
' general approach. 

The starting point is a microscopic, many body hamiltonian H written as the sum of a reference part Hq and an 
interaction term V. The interaction is assumed to be bilinear in some operator p(r), which is assumed either linear 
Q , in bosonic operators or quadratic in fermionic ones: 
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V = \ J d xdyp(x)w(x - y)p{y) (1) 

with a non singular (i.e. Fourier transformable) two body potential w. The properties of the reference system under 
the action of an external field h coupled to the order parameter p(r) are supposed known. No specific assumption on 
the reference system is made: in particular we do not need that Hq corresponds to a non interacting system, where 
Wick theorem applies (such a feature is crucial in setting up the QHRT equations). These requirements are indeed 
rather general and include several models of current interest in many body physics: quantum magnets (where p{r) 
represents the local spin variable), fermionic systems, like the Hubbard model, or even the Holstein model for the 
electron-phonon problem. 
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The first task is to build up a formal perturbative expansion of the partition function of the model: Z = 
Trexp(— f3H). Following a standard procedure 0, Z/Z Q can be written as the average over the reference distri- 
bution function of an imaginary time evolution operator U(/3). When this operator is written as a power series of the 
interaction w(r), we formally recover a perturbative expansion identical to that of a classical partition function for a 
(d+ 1) dimensional model. The additional "temporal" dimension is limited to the interval (0, 0) and w(r)5(t)//3 plays 
the role of classical two body interaction w c (r,t). The reference system of the associated classical model is implic- 
itly defined by requiring that its correlation functions coincide with those of the quantum reference hamiltonian Hq. 
Approximate mappings between a quantum model and an effective classical system have been proposed and studied 
in the literature || in order to clarify the role of thermal and quantum fluctuations. The novelty of our approach 
is that i) it is exact and ii) it applies to all temperature regimes, including the T — * limit. Having reduced the 
quantum problem to a classical one, we can directly apply the techniques developed in that framework. In particular, 
the already mentioned HRT is an implementation of the momentum space renormalization method which preserves 
information on the details of the microscopic hamiltonian. In HRT, different Fourier components of the two body 
interaction are included gradually, starting from the shortest wavelength: physically this corresponds to a smooth 
turning on of fluctuations over larger and larger lengthscales. This procedure can be carried out exactly by defining a 
sequence of auxiliary systems interacting via a potential whose Fourier components coincide with w(k) for k > Q and 
vanish elsewhere. As a result, we obtain a set of coupled differential equations expressing the change in the free energy 
and in the correlation functions of the model when a given Fourier component Q of the potential is included. The 
initial condition represents the system in which fluctuations are frozen, and in fact coincides with the known mean 
field result. As the cut-off wavevector Q — ► 0, the fully interacting system is recovered. Mean field approximation 
therefore corresponds to neglecting the change in the properties of the model as described by our differential equations. 
Details can be found in Ref. M . 

As an example, we study the the spin-S* antiferromagnetic Heisenberg model on a hypercubic lattice in d dimension: 

H = Hq + V = h e ig R 5 R + J S R- S R' ( 2 ) 

R <R,R'> 

where the sum is restricted to nearest neighbors and g is the antiferromagnetic wavevector of components gt = n. 
In applying the HRT approach to a quantum model, we decided to impose the cut-off Q only to the spatial Fourier 
components of w c (r,t). Different physical models might require other choices of the cut-off, whose only role is to 
continuously connect the reference system to the fully interacting one, which is recovered in the Q — > limit. The 
exact evolution equation for the Helmholtz free energy density a of the system describes how a is modified due to a 
change in the cut-off Q: 



+ ln[(l-fg(p,w)«;(p))(l+Fg(p',wMp))] } (3) 
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Here p' = p + g, w(p) = 2Jj(p) = 2J^ i cos(pi) is the Fourier transform of the interaction, the summation is over 
the Matsubara frequencies to n — 2im/P, the (d — 1) dimensional integral is restricted to the surface Sq defined by 
7(p) = — \J d 2 — Q 2 and the functions F®(p, =< S l (p, w)5 J '(— p, — uS) > are the Fourier transforms of the spin-spin 
dynamical correlation functions (in imaginary time) for a system where only fluctuations of wavevector p such that 
It (p*) I < \/ d 2 — Q 2 are included. The isotropy of the model implies F yy = F xx and F yx — —F xy . Analogous equations 
can be derived for the many spin correlation functions of the model. The (infinite) set of differential equations forms 
the QHRT hierarchy. When Q = d, fluctuations are neglected and the exact initial condition for the first QHRT 
equation (||) coincides with the mean field free energy density. The magnetic structure factors at Q = d can be 
explicitly written as: 

m z oj z + /Zj^ — w(p) z 
1 



m 2 uj 2 + /itj_ — w(p) 2 
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where m is the staggered magnetization and , are known functions of m. For S — 1/2: fj,± = 2Tm~~ 1 tanh~ 1 (2m) 
and = 4T(1 — Am 2 ) -1 . Note that the dependence of the transverse magnetic structure factors in (jj) on frequency 
and momentum is consistent with the first order spin wave result at zero temperature and reproduces the known single 
mode approximation which well represents antiferromagnetic correlations at low temperatures |^|,[l0| . The longitudinal 
correlations in equation (Q) are instead purely classical and, as such, satisfy the relationship T x\\{k) — Su(k) between 
longitudinal susceptibility Xll(^) an d the corresponding static structure factor Sn(k). 

Equation (|3|), although formally exact, is not closed because the evolution of the free energy depends on the unknown 
magnetic structure factors of the model F^(p, to). Therefore, we have to introduce some approximate parametrization 
of the structure factors in terms of the free energy. The simple approximation we have studied is to retain the form of 
Fjj(p, w) as given in Eq. (Q) but imposing thermodynamic consistency in order to determine the two scalar parameters 
fi± and fin. More precisely, for every Q we related the transverse and longitudinal staggered susceptibilities to the 
free energy via the exact sum rules: 

(/x_L - ZdJ)- 1 = Fg(g,uj = 0) = midaQ/dmy 1 

(Ml, - 2dJ)- x = Ff a {g,w = 0) - (d^Q/dm 2 )- 1 (5) 

From the adopted structure of the dynamical correlation functions, we also obtain the relationship between the 
parameters entering and the zero temperature non linear sigma model coupling constants: uniform transverse 



susceptibility xo — 1 / (4d) , spin wave velocity c = \J Adm and spin stiffness p s — m . The hydrodynamic relation 
Xoc 2 = p s is automatically satisfied by our ansatz for arbitrary spontaneous magnetization m. It is interesting to 
note that from our parametrization, the scaling of the spin wave velocity c on approaching the critical temperature 
(t = (T c — T)/T c — ► 0) gives c oc m that is c oc tr along the coexistence curve. The dynamic scaling hypothesis instead 
predicts c cx t v( - z ~^ where z is the dynamical critical exponent. By use of scaling laws and recalling that in our 
approximation the correlation critical exponent vanishes (77 = 0), we get z — d/2 which is the expected result for an 
antiferromagnet (i.e. model G) Equation (|J), together with (H|) and (||) give rise to a partial differential equation 
for the free energy density of the Heisenberg model a^(m) as a function of the cut-off Q and of the magnetization m. 
The frequency sum can be carried out analytically giving the final equation: 
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We numerically solved this partial differential equation for several values of the spin S and different temperatures in 
order to study the phase diagram of this system. Note that S just enters the theory through the initial condition 
a^(m) at Q = d, while the form of the differential equation is unaffected by S. Before showing the numerical 
results, however, it is useful to discuss the behavior of Eq. (||) near a phase transition. In particular, we studied 
the neighborhood of the critical point (in d > 2) and the low temperature region. Both at the critical point and 
along the coexistence curve the susceptibilities diverge due to the presence of critical fluctuations and Goldstone 
bosons, respectively. From Eq. (|^) we conclude that at long wavelengths (i.e. Q — > 0) and near a phase transition 
we have fx± ~ fi\\ ~ 2d J. In this region equation (^) simplifies and, by rescaling the free energy as a^Q~ d and the 
magnetization as mQ( 2 ~ d ' > / 2 1 it reduces to the RG equation obtained by Stanley et al. |l^] for a 0(3) symmetric 
</> hamiltonian. Such an equation has been analyzed near four dimension and proved to give the correct critical 
exponents to first order in the e = 4 — d expansion. In three dimensions, the numerical solution of the universal 
fixed point equation gives for the correlation length critical exponent the result v — 0.826, to be compared with the 
accepted value v — 0.71. The other critical exponents follow from the scaling laws, noting that our analytical form of 
the two point functions (|j) forces the anomalous dimension exponent to vanish r] = 0. We therefore find non classical 
exponents in three dimensions. Special care must be paid when dealing with the T — > limit of our equation. In 
this case, the asymptotic form of the equation changes and it can be shown to give rise, near a hypothetical quantum 
critical point, to critical exponents appropriate for a 0(3) model in d+ 1 dimensions as expected. A separate analysis 
should be carried out in the low temperature phase. If symmetry is spontaneously broken, following Chakravarty et al. 
p3[ , we may ask how quantum and thermal fluctuations modify the zero field magnetization. In order to answer this 
question we perform a Legendre transform on our equation (^): we first derive it with respect to the magnetization 
m obtaining an evolution equation for the magnetic field h^(m) at fixed m. Then we find the equation governing the 
evolution of the spontaneous magnetization rrfi implicitly defined by the requirement hP(rrfi) = at every Q. This 
procedure gives rise to a differential equation for rrfi . In the Q — > limit, taking into account that the longitudinal 
susceptibility diverges more slowly than Q~ 2 , QHRT reduces to a simple ordinary differential equation: 
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dmQ = K d f-^V 2 [tanh^mQ)]- 1 (7) 



where Kd is a geometrical factor (ratio between the solid angle and volume of the Brillouin zone). By introducing the 
rescaled variable g = \f4d{Q / \fd) d ~ 1 jvrfi , equation (Q) becomes identical, to order g 2 , to the known weak coupling 
RG equations for the non linear sigma model applied by Chakravarty et al. jl3| to the analysis of the antiferromagnetic 
Heisenberg model at long wavelengths and low temperatures. As an example, we plot in Fig. 1 the RG flux of g® 
obtained by the integration of the full QHRT equation (^). We clearly see the effect of the unstable zero temperature 
weak coupling fixed point while, for the nearest neighbor Heisenberg model, the other fixed point (g c ), governing the 
quantum critical regime, has no effect on the RG trajectories. This analysis shows that the single mode approximation 
to QHRT reproduces the correct long wavelength structure both near four and two dimensions. Furthermore, we expect 
QHRT to be superior to the weak coupling rcnormalization group equations because our non perturbative approach 
also describes the critical region and the high temperature regime where nrfi — > as Q — > 0, corresponding to g — > oo 
i.e. to the strong coupling phase of the non linear sigma model. 

Finally, we present few results of the numerical integration of Eq. (^) in three dimensions. As already pointed out 
in the classical case |(|, the HRT approach is able to correctly implement Maxwell construction at first order phase 
transitions and in fact the free energy density a^(m) at the end of the integration, i.e. in the Q — > limit, becomes 
rigorously flat in a finite region of the magnetization axis for T < T c . Therefore it is easy to extract from the numerical 
output the critical temperature and the coexistence curve, shown in Fig. 2 for several values of the spin S. The zero 
temperature limits of this curve agree within about 2% with the accepted estimates based on spin wave theory, Monte 
Carlo simulations or series expansions. Regrettably, for the spontaneous magnetization at finite temperature there 
are just few available results going beyond mean field approaches. Simulation data for the classical S* ^ oo case JlJ] 
and recent series expansion for the S = 1/2 model [ [l5[ seem to give somewhat larger coexistence regions. However, we 
believe that a more systematic analysis of these models by accurate numerical techniques is necessary before reaching 
a definite conclusion on the accuracy in the determination of the coexistence curve. The critical temperature for the 
classical model is known by several methods to be T c = 1.443 J while for the S — 1/2 case it has been recently 

estimated as T c = 0.946 J |J] by use of a newly developed Quantum Monte Carlo method and T c — 0.93 J Jll| by 
series expansions. Our results are a few percent lower, being T c = 1.419 J for S = oo and T c = 0.90 J for S = 1/2. 
From the solution of the QHRT equation we also obtain other important information on the model, for instance, 
the equation of state, the specific heat and also the temperature dependent dynamical structure factors, via analytic 
continuation of the adopted expressions (||). In order to improve the QHRT results we have just discussed, other 
approximate expressions for the magnetic structure factors should be examined, possibly keeping the same form (^) 
but allowing for a non trivial rcnormalization factor for the uniform susceptibility. This method can be applied in a 
straightforward way to other models of interest in quantum many body physics, like the Hubbard model, and may 
help to determine the location of the magnetic phase transitions and the possible occurrence of phase separation in a 
purely repulsive electron system. 
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FIG. 1. RG trajectories for the two dimensional Heisenberg model computed via numerical integration of the QHRT equation. 
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FIG. 2. Reduced spontaneous magnetization as a function of temperature for different values of the spin: S = 1/2 (triangles) 
5=1 (squares) S = 5/2 (circles) and S = oo (full line). 
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